Optimal transport by omni-potential flow and cosmological reconstruction 
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One of the simplest models used in studying the dynamics of large-scale structure in cosmology, 
known as the Zeldovich approximation, is equivalent to the three-dimensional inviscid Burgers equa- 
tion for potential flow. For smooth initial data and sufficiently short times it has the property that 
the mapping of the positions of fluid particles at any time ti to their positions at any time ti > ti 
is the gradient of a convex potential, a property we call omni-potentiality. Are there other flows 
with this property, that are not straightforward generalizations of Zeldovich flows? This is answered 
in the affirmative in both two and three dimensions. How general are such flows? Using a WKB 
technique we show that in two dimensions, for sufficiently short times, there are omni-potential 
flows with arbitrary smooth initial velocity. Mappings with a convex potential are known to be 
associated with the quadratic-cost optimal transport problem. This has important implications for 
the problem of reconstructing the dynamical history of the Universe from the knowledge of the 
present mass distribution. 

Dedicated to the memory of Roman Juszkiewicz 



I. INTRODUCTION 



Reconstruction in cosmology considers the following 
problem: one assumes that the present spatial distri- 
bution of masses (galaxies and clusters, including their 
dark-matter components) is known from observations, 
and one wants to reconstruct the dynamical history of 
the Universe all the way to the earliest epoch, when mat- 
ter and radiation decoupled (nearly 14 billion years ago). 
Peebles [22] introduced the reconstruction problem, and 
proposed a variational formulation for solving it on a rela- 
tively small spatial scale, that of the Local Group (which 
includes our own galaxy and neighboring ones). On scales 
much larger than that of the Local Group, which have 
been mapped in recent years through various projects 
such as the Sloan Digital Sky Survey [27], reconstruc- 
tion may be posed in the simplest cases as an instance 
of optimal mass transport. Indeed, Frisch et al. [13] 
showed that when the Zeldovich approximation [26] or 
a refinement thereof (cf. below) are applied to the rele- 
vant cosmological fluid equations, the correspondence be- 
tween the positions of mass elements initially (at decou- 
pling) and finally (at the present epoch) is the solution to 
an optimal mass transport problem with quadratic cost. 
This solution is uniquely prescribed by the marginals: 
the mass distribution at decoupling (essentially uniform) 
and its highly non-uniform present distribution. A strik- 
ing feature is that the sole knowledge of the current po- 



sitions of galaxies, without knowledge of their (proper) 
velocities, yields nevertheless a unique solution for this 
kind of large-scale reconstruction. 

It was then shown by Brenier et al. and by Loeper 
[7, 17] that, with prescribed marginals, unique recon- 
struction, not only of the Lagrangian map, but of the full 
dynamical history of matter elements, carries over to the 
Euler-Poisson model, whose validity extends much be- 
yond that of the Zeldovich approximation. Its unique so- 
lution is again obtained from an optimal transport prob- 
lem with a convex cost function, expressible as a space- 
time integral of a suitable action, a problem whose nu- 
merical resolution remains a challenge. 

As is well known, the mass transport problem was in- 
troduced by Monge [19] more than two hundred years 
ago, and the theory took its modern shape after the 1942 
work of Kantorovich [16] (see, e.g., Villani [25] for re- 
view). 

The Zeldovich approximation [26] was introduced in 
1970 as a first formulation in terms of Lagrangian coordi- 
nates of the growth of density perturbations. It replaces 
the full Euler-Poisson equations by basically the three- 
dimensional inviscid Burgers equation (written here in 
standard fluid dynamical notation) 
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The validity of the Zeldovich approximation is con- 
trolled by how close one is to decoupling, but in a scale- 
dependent way: at very large scales, the Zeldovich ap- 
proximation remains valid up to the present epoch; at 
very small scales, the formation of multi-stream caustics 
quickly ruins not only the validity of the Zeldovich ap- 
proximation, but even that of the Euler-Poisson model. 
An immediate consequence of (1) is that the velocity of 
any fluid particle remains constant in the course of time 



2 



and that the trajectories are straight lines. We denote by 
q the initial (Lagrangian) fluid particle positions and by 
x their (Eulerian) positions at the current epoch t = T. 
The Lagrangian map associated with the Zeldovich ap- 
proximation is 

q^x = V q (^-+TM<l)) , (2) 

where <po{q) is the initial velocity potential. For suffi- 
ciently small T and a sufficiently smooth initial poten- 
tial, the Lagrangian map is thus the gradient of a con- 
vex function, a property shared by the next-order ap- 
proximation, which will be discussed in Sec. V. This 
is why reconstruction is linked to optimal transport; in- 
deed, a theorem of Brenier [5] states that the solution 
to the Monge optimal transport problem with quadratic 
cost is a gradient of a convex function, which satisfies 
a Monge- Ampere equation. The method of cosmological 
reconstruction in which one assumes that the Lagrangian 
map has a convex potential and then numerically solves a 
quadratic-cost optimal transport problem (after suitable 
discretization) is called the Monge- Ampere-Kantorovich 
(MAK) method [13]. 

The Zeldovich approximation gives us some in- 
sight into the full temporal history of mass elements. 
An important consequence of (1) is that for any 
< t\ < t 2 < T the mapping of particle positions at time 
t\ to their positions at time t 2 is also a gradient of a con- 
vex function. When the flow-induced mapping between 
any two times is potential, the flow is here called omni- 
potential. As we shall see, the velocity field associated 
with such flows has the property of being simultaneously 
potential in Eulerian coordinates (in cosmology, this con- 
straint stems from the potential character of gravity and 
the expansion of the Universe), as well as in Lagrangian 
coordinates (which allows reconstruction by solving an 
optimal transport problem). 

We are of course led to ask whether there exist omni- 
potential flows other than Zeldovich/Burgers ones or 
trivial variants thereof. Investigating this issue is the 
central topic of our paper. 

In Section II we show that omni-potentiality can be 
rccxpressed geometrically and algebraically in terms of 
Hessian matrices and recast as a set of one or several 
partial differential equations (depending on the space di- 
mension). In Section III we use an algebraic method to 
construct explicit non-trivial examples of omni-potential 
flows in two and three dimensions. These are rather spe- 
cial and we are led to ask how general are omni-potential 
flows. In Section IV we construct a fairly general class of 
two-dimensional omni-potential flows, leaving the prob- 
lem open in higher dimensions. In Section V we re- 
turn to questions of cosmological relevance: To what ex- 
tent are the full solutions to the cosmological equations 
omni-potential? Why is the validity of MAK reconstruc- 
tion better than that of the Zeldovich approximation, as 
pointed out by Mohayacc ct al. [18]? In Section VI we list 
some open problems and make concluding remarks. Fi- 



nally, in the Appendix we characterize sets of commuting 
symmetric matrices by constructing suitable invariants. 

II. CRITERIA FOR OMNI-POTENTIALITY OF 
FLOWS 

In the present paper we study the kinematics of omni- 
potential flows. We start by recalling the basic definitions 
of the Lagrangian and Eulerian description of a flow — 
a motion of fluid regarded as a continuum of infinitesimal 
fluid particles (whose mathematical abstraction is a point 
particle). 

Denote by v(x,t) the velocity of the fluid measured 
at point x in space and at time t. It is usually called 
the Eulerian velocity (the velocity measured at a fixed 
position in the laboratory frame). The motion of a fluid 
particle satisfies the ordinary differential equation 

x = v(x, t), 

which has to be supplemented by the initial condition 
as|t=o = <?• 

If the velocity field v(x,t) is prescribed and sufficiently 
regular, one can solve this initial value problem, at least 
locally in time, and obtain the mapping q i-> x(q,t) 
called the Lagrangian map. It takes a particle at the 
Lagrangian position q and carries it to the Eulerian po- 
sition x. For a fixed q, the curve x(q,t), parameterized 
by varying time t, is the trajectory of the particle, whose 
Lagrangian position is q. When, in a field associated 
with the flow, we perform the substitution x — > x(q,t), 
we obtain its Lagrangian description, in which the field 
is now a function of the Lagrangian coordinates q. For 
instance, v(x(q, t),t) is called the Lagrangian velocity. 

In this paper, we consider flows v(x,i) defined, for 
simplicity, in the entire space M. d , but restricted to a finite 
time interval [0,T]. The flow induces a set of mappings 
of space: given two arbitrary times t and r, such that 
< t < t < T, the mapping from fluid particle positions 
at time t to their positions at time r is here called the 
(t, -remapping. The (0, remapping is just the standard 
Lagrangian map. 

As stated in the Introduction, for any two times t and 
r, such that < t < r < T, the (t, r)-mappings induced 
by omni-potential flows are required to be the gradients 
of a convex potential $(q, t; r): 

q^x = V,$(9,t;T). (3) 

Such a mapping is here called potential and 4>(<7, t;r) is 
called the (t, r)-potential. Given any three times to, t 
and r such that 0<t n <t<T<T, the (to, i)-mapping 
composed with the (t, remapping obviously yields the 
(to, remapping. This semigroup associativity, combined 
with omni-potentiality, implies 

$( q ,to;T) = $(V,$(g,t ;t),t;T), (4) 
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FIG. 1: Potential composition of two potential maps. 

which is illustrated in Fig. 1. 

Differentiation of (3) with respect to time r shows that 
the velocity, v(q, t;r), is also potential: 

v(q,t;r) = V q (p(q,t;r), 

where 

ip{q,t;r) = d T <$>(q,t;r). 

In this notation, v(q, 0; t) = V q (p(q, 0; t) is the La- 
grangian velocity of a fluid particle, while its Eulerian ve- 
locity is v(x,t;t) — V x f(x,t;t); thus, in omni-potential 
flow, the velocity is potential in both Eulerian and La- 
grangian coordinates. In three dimensions, this is equiv- 
alent, locally, to the statement that the vorticity (the 
curl of the velocity) vanishes in both coordinate systems. 
We shall call this dual potentiality in both Eulerian and 
Lagrangian coordinates bi-potentiality. 

We demand that all (t, r)-potentials of omni-potential 
flows be convex in the space variables. This implies that 
the associated (t, r)-mappings have inverses that are also 
potential, the potentials associated with the inverse map- 
pings being the Legendre transforms of those of the di- 
rect mappings (see Ref. [7], Appendix C and references 
therein). Actually, for a sufficiently smooth flow, the con- 
verse is true: invertibility implies convexity of the maps, 
as now explained. Let us denote by H(f) the Hessian 
of a (twice differentiablc) function f(q), i.e., the d x d 
matrix with the entries 

««(/) '>;,„!■ 

We now assume that the potentials are twice differen- 
tiable in the space variables and that their Hessians 
H($(q, t; t)) depend continuously on their time argu- 
ments. We then observe that, for coinciding times, the 
(t, t)-mapping is the identity mapping, clearly having the 
convex potential 

*{q,t;t) = \q\ 2 /2, 



whose Hessian is the identity. As the two times separate, 
loss of convexity would require one or several eigenval- 
ues of the Hessian to change sign and thus to go through 
zero; at such an instant, the Jacobian matrix (which for 
a potential mapping coincides with the Hessian of the 
potential) becomes degenerate; then, generically, the in- 
verse mapping ceases to exist, i.e., the property of the 
(t, -remapping to be bijective gets lost (in the cosmo- 
logical context this amounts to shell crossing leading to 
-multi- streaming) . 

Below, we establish criteria for omni-potentiality. In 
Sec. II A we prove that a flow is omni-potential, when- 
ever the Hessians of the potentials $(q, t; r), calculated at 
any two points of a trajectory, commute. In Sec. II B we 
present an equivalent condition for the commutativity of 
Hessians: along any trajectory, at any time t, the Hessian 
TL(t) = H($>(q;0,t)) and its time derivative should com- 
mute; this is used to show that omni-potentiality is equiv- 
alent to potentiality of both the Lagrangian and the Eu- 
lerian flow velocities (plus the convexity constraints). In 
Sec. II C we discuss simple examples of omni-potentiality: 
Zcldovich and Zeldovich-type flows. Finally, in Sec. II D 
we derive a partial differential equation for the poten- 
tial of a two-dimensional omni-potential flow. It states 
that a suitable expression constructed from second-order 
derivatives depends only on the Lagrangian coordinates, 
but not on time. In other words, the Hessian of the poten- 
tial possesses an invariant, whose value depends only on 
the trajectory. We construct similar invariants in higher- 
dimensional spaces in the Appendix. 



A. Commutation of Hessians of the potential 

The semigroup associativity (4) involves the compo- 
sition of two potential maps. We shall show that, in 
general, such a composition is potential if and only if 
the Hessians commute. Basically this stems from the 
well-known theorem that the product of two symmetric 
matrices is symmetric only when they commute. The 
problem we are now addressing is illustrated in Fig. 1, 
which sketches the action of three mappings along the 
same trajectory. We observe that the (t, remapping is 
the composition of the inverse of the (to, t)-mapping with 
the (to, -remapping. We shall now show that its poten- 
tiality is equivalent to the commutation of the Hessians 
of the (to, t)-mapping and of the (to, remapping. 

Specifically, we assume that, for any times t and r, 
such that < t < r < T, the (t, r)-mapping (3) of R d 
into itself is a bijection that, together with its inverse, is 
smooth (i.e., has as many continuous derivatives, as we 
might need). We denote 

$i(g) = $(q, t ; t), $ 2 (<?) = $(<?, t), 

where $(q,t ;t) and ^(q,t ;r) are the potentials of 
the (to, t)-mapping and the (to, remapping, respectively. 
The required potentiality of the (t, -remapping implies 
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that the Jacobian matrix 



is symmetric: 



for any pair of indices 1 < i, j < d. The converse is also 
true, at least locally in space. By the chain rule, 

H mn (^ 2 ) = d q J m 



fe=l 
d 

$^£m«fen(*l)- 



k=l 



Therefore, 



da 



dx 



^($ 2 )H" 1 ($i), 



(5) 



where denotes the matrix inverse to For 

the r.h.s. of (5) to be a symmetric matrix, the matri- 
ces H{<&2) and must commute, which is equiva- 
lent to the commutation of the two Hessians ^($2) and 

Reversing the arguments, we establish that commuta- 
tion of Hessians $(q, to; t) along a trajectory, 

H($(q, t ; t))U(<5>(q, to; r)) = H($(q, t ; r))H(^(q, t ; t)), 

at any times t and r is necessary for the mapping to 
be potential. In particular, commutation of Hessians of 
the potential $(<7, 0; t) of the Lagrangian map along each 
trajectory (for each fixed q), together with invertibility 
or convexity, is equivalent to omni-potentiality. 

By the theorem on codiagonalizability of commuting 
(real) symmetric matrices (see, e.g., Ref. [15], pp. 50-51) 
the equivalent condition is that the Hessians of the po- 
tential $>(q,to;t) calculated at different times t for the 
same coordinate q and the same to are codiagonalizable, 
i.e., can be transformed into the diagonal form using the 
same unitary matrix. In other words, along any trajec- 
tory, only the eigenvalues of the Hessian of the potential 
but not the eigendirections are allowed to vary. 

So far, we have shown the commutation — along a 
given trajectory — of the Hessians of the potentials of 
the (£o,i)-mapping and the (to, remapping for the same 
starting time to (e.g., for the (to, To)-mapping and the 
(to, ti)-mapping of Fig. 2). A similar argument proves 
the commutation of the Hessians of the potentials of two 
mappings, such that the ending time of one of them co- 
incides with the starting time of the second one (e.g., 
for (to,ti)-mapping and the (t\, ri)-mapping). Combin- 
ing these two results and relying again on the theorem 
on codiagonalizability of symmetric commuting matrices, 
we find that the Hessians of the potentials of the (to, to)- 
mapping and of the (ti, ri)-mapping commute. Thus we 
have established that, along any given trajectory, the 
Hessians associated with arbitrary pairs of times com- 
mute. 




FIG. 2: Commutation of Hessians of the potential along a tra- 
jectory. A sketch of a trajectory and flow-induced mappings 
from times to and ti to times to and n . 



B. Commutation of Hessians and their time 
derivatives: bi-potential velocities 

Here we give an alternative formulation of omni- 
potentiality in terms of commutation of Hessians and 
their time derivatives. We need some preparatory ma- 
terial regarding d x d symmetric matrices with smooth 
time dependence. Suppose at any two times t and t' they 
commute: 



H(t)H(t') = H(t')H(t). 



(6) 



Differentiating this equation in t' and letting t' = t, we 
find that at any time H (t) commutes with its time deriva- 
tive H(t): 



H(t)H(t) = H(t)H(t). 



(7) 



We shall show now that the converse is also true: if at any 
time t (i) relation (7) is satisfied, and (ii) all eigenvalues 
Aj of the symmetric matrix H(t) are distinct, then (6) 
holds true for any times t and t\. 

Since H(t) is symmetric, it can be expressed as 



H(t) = U*(t)A(t)U(t), 



(8) 



where U is a unitary d x d matrix, and A is diagonal. 
Consider the identity U (t)U l (t) — I, where I is the iden- 
tity matrix; differentiating it in time yields 



UU 



7* = -UU* = - (UU^ 



Thus X = UU 1 is an antisymmetric matrix. 

We substitute (8) into (7), and multiply it by U on the 
left and by £/' on the right, obtaining 

AAA - A 2 A = AA 2 - AAA, 

i.e., the matrix AA — AA commutes with A. This ma- 
trix is symmetric due to the antisymmetry of A. By the 
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theorem on codiagonalizability of commuting symmetric 
matrices, the matrices XA — AX and A are thus simul- 
taneously diagonalizable; but A is already diagonal, and 
hence so is XA — AX. The entries of the latter matrix are 
(Xj — Xi)Xij, and therefore Xij = for all i ^ j (we use 
here the condition that all eigenvalues of H are distinct). 
The antisymmetry of X implies that all diagonal entries 
of X also vanish, and thus X = 0. Therefore, 

U = -XV 

vanishes. In other words, variation of H in time consists 
solely of variation of its eigenvalues Aj. This implies (6). 

The restriction that all the eigenvalues be distinct does 
not significantly affect the generality of the statement: if 
they are not distinct at some isolated times, the relation 
U = remains satisfied at these times by continuity. 

Returning to the problem of omni-potcntiality, we now 
take 

H{t) = «($(qf,0;t)) = H. 

This Hessian H is the Jacobian of the Lagrangian map 
V$(q, 0;t). By the above result, omni-potentiality is 
equivalent to the commutation, at any time, of the Jaco- 
bian and its time derivative. 

We have seen earlier that omni-potential flow has a 
velocity which is bi-potential, i.e., potential in both Eu- 
lerian and Lagrangian coordinates. The statement just 
derived allows to prove the converse, namely that a flow 
with a bi-potential velocity (and some convexity require- 
ments) is omni-potential. Let us denote by v L (q, t) and 
v E (a:,i) the Lagrangian and Eulerian velocity, respec- 
tively. Since x{q,i) is the Lagrangian map, we obviously 
have 

v E (x,t)=^(q(x,t),t), (9) 

where q(x,t) is the inverse Lagrangian map, whose Ja- 
cobian is W.^ 1 . We now calculate the Eulerian velocity 
gradient, using (9). By the chain rule, for any i and j, 
we have: 

d 

d Xi vf{x,t) = 9 qm vf(q, t) 

m—l 

d 

m— 1 
d 

^ ^ )im%mj- 
m—l 

Thus, the Eulerian velocity gradient is the product of 
the matrices Hr 1 and H. For the Eulerian velocity to be 
potential, it is necessary and (locally) sufficient that this 
product be a symmetric matrix. The commutation of the 
symmetric matrices V.^ 1 and % is equivalent to the com- 
mutation of H and H. Equivalence to omni-potentiality 
follows from the statement derived above. 



C. Zeldovich and Zeldovich-type flows 

In the Zeldovich approximation each particle keeps its 
initial velocity unaltered in the course of time, and hence 
particles move along straight lines (at least before multi- 
streaming occurs). The Lagrangian map at time t is 

q^x = V q (\?L+ttpo(q)j , 

where <fio(q) is the velocity potential, prescribed at t — 0. 
The Hessian of this map is I + tH(ifo), where I is the 
identity matrix and the matrix H(ipo) is the Hessian of 
the initial potential. For a given q, the eigendirections 
of the associated Hessian are those of the Hessian of the 
initial velocity potential. Clearly, all these Hessians com- 
mute and, by the results of Sec. II A, such a flow is omni- 
potential. 

More general examples of omni-potential flows can be 
constructed by performing an arbitrary nonlinear trans- 
formation of the time and by time-dependent zooming 
of space. In space of any dimension d > 2, consider the 
flows defined by the potentials 

*(9,0;t) = M*)^+»K*to>(<l), (io) 

where fi(t) and r](t) are arbitrary functions of time. 
Clearly, these are again omni-potential. 

In general, the trajectories associated with (10) are 
not straight lines. However, if we look at them with a 
time-dependent magnifying glass which applies a zoom- 
ing factor l/fi(t), they become straight. Furthermore, if 
we introduce a new time variable t' = rj{t)/ ' [i{t), particles 
move again with a constant velocity. Hence, the flows de- 
fined by (10) are trivial generalizations of Zeldovich flows, 
and will here be called Zeldovich-type flows. 

Our goal is to find omni-potential flows that are not of 
this type. 



D. A linear second-order PDE for two-dimensional 
omni-potential flow 

We derive here a differential equation for the potential 
of a two-dimensional omni-potential flow. It turns out to 
be a linear second-order PDE. 

Consider a symmetric 2x2 matrix H. Suppose its 
eigenvector associated with eigenvalue A makes angle 6 
with the cartesian axis q\. Thus, 

Hu cos 9 + H 12 sin 9 = A cos 9, 

H12 cos 9 + H 2 2 sin 9 = A sin 9. 

In order to eliminate the eigenvalue A, we multiply the 
first of these equations by sin 9 and the second one by 
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cos 9. Subtracting afterwards the second equation from 
the first one, we obtain 



Hu — H 2 



H 



cot 29. 



12 



(11) 



Prescribing the r.h.s. of (11) uniquely defines the orthog- 
onal frame of the two eigendirections. (The values of 
cot 29 define the angle 9 modulo 7r/2; however, chang- 
ing 9 — > 9 + it/ 2 swaps the eigendirections, but does not 
affect the set of eigendirections.) 

In an omni-potential flow, the eigendirections of the 
Hessians of the (0, t)-potentials should depend only on 
the Lagrangian position q and not on the time t. Let 
<&(g, t) be a two-dimensional omni-potential flow and let 
us denote by g{q) the common value of cot 29 along the 
particle trajectory emanating from q. It then follows 
from (11) that 



(d 2 -d 2 )<I> = q(q)d 2 



(12) 



The search for two-dimensional omni-potential flow 
has thus been reduced to finding solutions to (12) for 
suitably prescribed functions g(q). 



III. EXAMPLES OF OMNI-POTENTIAL 
FLOWS IN TWO AND THREE DIMENSIONS 

The main question that we address in this paper is 
whether omni-potential flows exist that are not of Zel- 
dovich type. In this section we give a positive reply to this 
question both in the two- and three-dimensional spaces 
by providing explicit examples of polynomial potentials 
for mappings induced by such flows. 

An example of an omni-potential flow in a space of ar- 
bitrary dimension is provided by spherically-symmetric 
potentials of the form $(|q|,t): a simple calculation re- 
veals the commutation of Hessians of such potentials, 
calculated at different times at different points of a tra- 
jectory. This example shares with Zeldovich flows the 
property that the trajectories are straight lines — in this 
case, in the radial direction. We would like to construct 
less symmetric examples. 



of \q\ 2 /2 (which is always a solution to (12)), then the 
flow with the potential \q\ 2 /2 + ai(t)$i(q) + a 2 (t)$ 2 (q) 
is omni-potential, and is typically not of Zeldovich type; 
for this, the functions of time ct\ and a 2 must be linearly 
independent and sufficiently small, so as not to spoil the 
convexity stemming from the \q\ 2 /2 term. 

When g(q) is a ratio of homogeneous polynomials of q 
(say of the same degree to), solutions to (12) can be ob- 
tained by a purely algebraic method. A solution can be 

(2) 

sought in the form of a homogeneous polynomial, p n (q), 
of degree n > m + 2; then (12) reduces to a system of 

TO+n— 1 equations for the coefficients of p^ (<?) and g(q). 

(In what follows, plf* denotes certain homogeneous poly- 
nomials of degree n defined in R d .) The function g(q) 
involves 2to + 1 independent coefficients (since the nu- 
merator and denominator can be multiplied by any con- 

(2) 

stant without changing g(qj), and the polynomial p n (q) 
involves n independent coefficients (since a solution to 
(12) can be multiplied by any constant without yielding 
a new independent solution) . Comparison of the number 
of equations, to + n — 1, with the total number of the 
unknown coefficients, 2to + n + 1, suggests that we can 
construct a family of such solutions, parameterized by 
to + 2 coefficients of g(q). However, the system of equa- 
tions for the coefficients is, in general, nonlinear, and 
hence its solvability cannot be established just by count- 
ing the numbers of the unknowns and equations. When 
g(q) is the ratio of linear functions, the equations for the 

coefficients of Pn\q) are linear, and can be solved for 
any prescribed coefficients of g(q). 

Since the potential $ is required to be a convex func- 
tion on the entire plane K 2 , we start by seeking homo- 
geneous polynomials Pn\q) involving only even powers 
of qi and q 2 . An instance of a solvable linear system of 
equations yielding the coefficients of such polynomials is 
obtained for 



qgi - bq 2 

Q1Q2 



(13) 



where the coefficients a and 6 may take arbitrary pre- 
set values. For such g(q), a homogeneous polynomial of 
degree 2k, k > 2, solving (12), is 



A. Particular examples of two-dimensional 
omni-potential flow 

In M 2 , the problem of finding omni-potential flows has 
been recast into the form of the partial differential equa- 
tion (12) with the initial condition \q\ 2 /2 (which gener- 
ates the identity map) . We can therefore construct exam- 
ples of two-dimensional omni-potential flows by finding 
different solutions to (12) for a given function g(q) in the 
r.h.s. By linearity, any linear combination of such solu- 
tions with time-dependent coefficients is also a solution to 
(12). For example, if $i(q) and $2(9) are two sufficiently 
smooth independent solutions that are also independent 



k I i-1 

P2IV92) = £ l[(2k - 1 + 2j(a - 1)) 
*=o \j=o 



k-l-i 

x H (2fc-l + 2j(6-l)) 

3=0 



k\ qfql 



2i 2(fe-i) 



i\(k-i)\(2k- 1) 



• (14) 



In particular, the first low-degree polynomial solutions 
are: 

pi 2 V<? 2 ) = (2a+l)q 4 1 +6q 2 1 q 2 2 + (2b+l)q 4 2 , (15) 

?4 2) (91,^2) = (4 a +l)(2a + 3)^ + 15(2a + 3) gig2 2 

+15(26 + 3)q 2 ql + (46 + 1)(26 + 3)^. (16) 
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As one can see, the polynomial (14) vanishes identically 
for integer j > 1 and j > 1 such that j + j < k — 1 , and 

2fc-l - a 2fe — 1 



2= 1 



2j 



6=1- 



2j 



(17) 



For these isolated values in the plane of parameters, there 
exist nevertheless two independent solutions, namely 



d_ 

da 



(2) 
P2k 



and 



a— a. b—b 



d_ 

db 



(2) 
P2k 



a— a, b—b 



This can be easily seen by differentiating (12) in a and b 
and substituting the parameter values (17). 

Clearly, (q) is convex, if all coefficients in (14) arc 
positive, i.e., if 

min(a,6) > -l/(2fc - 2). (18) 
Thus, the potentials 

$(<?, t) = » 2 (t) J^- + £ /x 2fc (t)p< 2 fc > ( 9l , g 2 ) (19) 

fc>2 

are convex for min(o, b) > 0, if in addition all /i2k(t) are 
non-negative (this condition is sufficient, but not neces- 
sary) and tend to zero fast enough to guarantee conver- 
gence of the series at any point q and termwise differ- 
entiability of (19) in the spatial variables. If the sum 
(19) is finite and the maximum degree of the polyno- 
mials involved is 2K, then the potential is convex for 
min(a, b) > —1/(2K — 2). The initial condition is satis- 
fied provided \i 2 (0) = 1 and [i 2 k (0) = for all k > 1 . The 
convex potentials (19) satisfy all requirements for omni- 
potcntial flows in the plane, and are not of Zcldovich 
type. 

So far, we have considered only even-degree homo- 
geneous polynomial solutions or linear combinations 
thereof. Is an admixture of odd-degree homogeneous 
polynomials permitted? If such an odd-degrce addition 
has a degree higher than that of the highest even-degree 
homogenous polynomial comprised in the solution, then 
convexity in the whole plane is ruled out. However, in a 
finite domain, convexity need not be lost, provided the 
odd polynomial is scaled by a sufficiently small factor. 
This is precisely what happens when g(q) is given by 
(13): a homogeneous polynomial of odd degree 2k + 1 
can be a solution to (12) only for k > 2 and 

a = b= -l/(fc-l). (20) 

The solution for these parameter values is 



P2fe+i(9i,92) = C!p k {qi,q2) + c 2 pk(q2,qi), 



where 



fe-i 

t»fc(«i,ffi»)=5z ( n 

i=0 \j=0 



(2(fc-j) + l)(fc-l- j ) 
(j + l)(2j-l) 



2i 2(fe-i)+l 

x 91 9 2 



and ci and c 2 are arbitrary constants. For instance, for 
a = b = —1 fifth-degree polynomial solutions are 



c 2 q\- 



Comparison of conditions (18) and (20) shows that the 
degree of an odd-degree polynomial solution for g(q) de- 
fined by (13) is higher than the degree of any convex 
even-degree polynomials existing for the chosen g{q). 



B. Examples of omni-potential flow in dimension 

d > 3 

The approach that has been applied in the previ- 
ous subsection for construction of an example in dimen- 
sion two cannot be immediately generalized to higher- 
dimensional spaces: while in R 2 a single invariant deter- 
mines whether two symmetric matrices are codiagonal- 
izable and hence equation (11) fixes the set of eigendi- 
rections of a symmetric matrix H, in M 3 at least three 
such invariants must be considered simultaneously (see 
the Appendix). In dimension three, equations (A7)-(A9) 
applied for the entries of the Hessian of an unknown po- 
tential give rise to three partial differential equations in 
the potential, 



d 2 $ / (d 2 - d 2 )$ d 2 $ d 2 

91)93 I ^9l)9l 92,92/ , 92,93 91,93 



d 2 $ 

"9l,92^ 



d 2 $ 

"9l ,92 ^ 



+ 



d 2 $ d 2 $ 

"91)93 "92,93 . 



/ (d 2 -d 2 )$' 

= Ufa) + ■ ^ 



d 2 $ 

"91,92^ 



(^92,92 -^3,93) $ + 9 9 2 l,93 $ 

(d, 



d 2 $ 

"92,93 



d 2 $ d 2 $ 

91,92 "91,93 , 



(21) 



d 2 $ / 

92, 93 / / 

92-^ \9dQ) 

"91,93 \ 



52 ^^2 ^ 

91 92,92/ 

<9 2 $ 

"91,92^ 



(a 2 - s 2 )$ 
= . 92(g) -,g 3 (g)^ q ^ L - 



d 2 $ 

"91,92^ 



(22) 



.93(9) 



' (d 2 -d 2 )$ <9 2 $ a 2 

V 92,92 93,93 , 91,93 91,92 



a 2 $ 

"92 ,93 



+ 



o 2 $ a 2 $ 

91,92^ "91,93 , 



(a 2 - a 2 

_ V 93 ,93 



<9 2 $ 

91,93 



,$ a 2 $ 9 2 $ 

91 ,9l/ _|_ 91,92 92 ,93 



<9 2 $ d 2 

92,93 91,92 



(23) 



that must be satisfied simultaneously. Here the time- 
independent quantities gt (q) are related to the invariants 
7 2 i' fc ^ introduced in the Appendix: 

3i(9) 



(3,1) (3,3) 
721 + 721 , 



I \ (3,2) , , 

92(q) = 7 2 i + !) 



93(g) 



(3,3) 
721 • 



The invariants are combinations of the entries of a sym- 
metric matrix, which depend only on the set of eigendi- 
rections and not on the eigenvalues. The algebraic the- 
ory of the invariants presented in the Appendix does not 
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take into account the properties of the Hessian, stemming 
from the specific structure of its entries (for instance, 
each row and a column of the Hessian is a gradient, which 
implies certain differential relations between the entries). 
It is unclear how to prescribe gk(q), taking into account 
these additional constraints, for the equations to have at 
least two distinct solutions. 

Because of this difficulty, instead of considering the in- 
variants and solving equations (21)-(23), we exploit the 
fact that omni-potentiality of flows amounts to commuta- 
tion of the Hessians of the various (ti, i 2 )-mappings along 
any trajectory (see Sec. II A). We shall construct our ex- 
amples in R d using the following strategy. The potential 
is sought in the form of a linear combination of "build- 
ing blocks" with time-dependent coefficients. One of the 
building blocks is prescribed; we take it to be a homo- 
geneous polynomial, Pm\q), of degree to. All the other 
building blocks are then required to have their Hessians 
commuting with that of the prescribed building block. 
The function \q\ 2 , whose Hessian is the identity matrix, 
constitutes a trivial solution. We can try finding other 
building blocks in the form of homogeneous polynomials 

Pn\q) °f some degree n. If we succeed, it is easy to show 
that any linear combination of such building blocks (with 
the convexity restriction) will define an omni-potential 
flow. We seek such polynomials by requiring the vanish- 
ing of the non-diagonal entries of the commutator of the 
two Hessians, viz. 

CfoWpW) = H{p<$)H{pW) H(pW)H(p%>). (24) 

Unfortunately, in general this strategy does not work, 
as now explained. A homogeneous polynomial of degree 
n has 

{n + d- 1)! 
n\(d- 1)! 

coefficients. Non-diagonal entries of the commutator C 
are homogeneous polynomials of degree to + n — 4. The 
commutator is antisymmetric (recall that the Hessians 
are symmetric matrices), hence we have to consider the 
d(d — l)/2 non-diagonal entries of C. Thus, in general, 
we have to solve 

d(m + n + d — 5)! 
2{m + n - 4)\(d - 2)\ 

equations, a number which is easily seen to exceed the 
number of coefficients, 

(m + d-l)\ {n + d-l)\ 
m\(d- 1)! + ' 

So, the problem is overdetermined. 

Nevertheless, potentials having all the required prop- 
erties can be constructed in R d (d > 2), if all the build- 
ing blocks are restricted to be homogeneous polynomi- 
als symmetric in all their arguments (i.e., invariant un- 
der any permutation of the spatial variables % «-> qj). 



Such building blocks have the following significant ad- 
vantage: it suffices to consider any of the polynomial 
equations arising from non-diagonal entries of the com- 
mutator (24) (referred to as "commutator equations") - 
all these equations are equivalent by virtue of the sym- 
metry. In what follows, we implement this "symmetric 
building block strategy" in two cases, in R d for d > 3 
with just one unknown building block, and in R 3 with 
infinitely many ones. 

Now, we focus on the symmetric polynomials 

A d \q) = E^ 4 + ?EE^ 2 , (25) 

= E^+«EE^ 2 

i—1 i—1 j—1 

+ b E <W (26) 

l<i<j<k<d 

(When d = 3, the last sum in p$\q) reduces to a single 
term bq^q^q 2 ) We consider polynomials involving only 
even powers of the spatial variables qj, because we seek 
solutions that are convex functions. 

To be specific, we consider the commutator equation 
Ci 2 {p { i\pt ] ) = in R d , d > 3. The l.h.s. is a polynomial 
of degree 6. Since in p^ and p^ any power of qi and 
<72 is even, C\i is proportional to q\q 2 , and every variable 
enters into the polynomial Cyij '(qiqi) only in an even 
power. Since both potentials are symmetric in qi and 
<Z2, C12 = for q\ = q 2l and hence C\i is divisible by 
q\ — q\. The polynomial Ci2/(<7i32(? 2 ~ q 2 )) i s 01 the 
second degree; it is thus just a sum of q? with certain 
coefficients. Because the potentials are symmetric in qj, 
it has the form 

d 

a i{ql + ql) + a 2E < ?j- 

Hence, we have three independent parameters, a, b and 
c, and two equations to satisfy. Calculating the coeffi- 
cients «i and «2 and letting them vanish, we find that 
the Hessians of p^ (q) and p^ (q) commute for 

a = 15c/(12-c), (27) 
b = 75c 2 /((I2-c)(3 + c)). (28) 

Thus, the potential 

$( 9 ,t) = /i 2 (t)^ +^{t) P ( t\q) +^{t)pf\q) (29) 

defines a non-Zcldovich-type omni-potential flow in R d 
(d > 3). Polynomials pf\q) and p^\q) are convex pro- 
vided < c < 12; hence, if all Hi(t) > 0, potential (29) 
is convex for c from this interval. For c ^ 2, P^\q) and 
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Pq do not possess spherical symmetry, and hence the po- 
tential (29) is not spherically symmetric. Restrictions of 

pf\q) and p^\q) onto the plane q 3 = ... = q d = coin- 

(2) (2) 

cide with the polynomials p\ (qi,q 2 ) (15) and p 6 (qi , q 2 ) 
(16)fora = &=(6-£)/(22). 

Henceforth, for the sake of simplicity, we assume that 
the problem is three-dimensional. In the remainder of the 
subsection we shall implement the building block strat- 
egy with unknown blocks that are arbitrary even-degree 

polynomials p^nil) f° r n > 2. The polynomial p^\q) 
(25) for d = 3 remains our prescribed building block. 
By the theorem on codiagonalizability of symmetric ma- 
trices, commutation with the Hessian of p 4 (q) implies, 
that the Hessians of any two polynomials from this fam- 
ily commute. (This is taking place generically, i.e., at 
those points in R 3 , where the Hessian of p^ (q) does not 
possess equal eigenvalues; at non-generic points the com- 
mutation follows from continuity of the Hessians and the 
commutation at generic points, which are present in any 
neighborhood of a non-generic point.) 
The polynomial 

P ( il{Q)= E ^cifq^qf (30) 

i+j+k=n 

is symmetric, whenever 

Oi j k does not depend on the order of subscripts i, j, k. 

(31) 

Straightforward algebra yields 

Ci2(p&U 3) ) = 8 gi <6 E ",.,.,'/?' \V \V 

z,i,fe>o 

i-\-j-\-k=n 

x (y(c-6)(g?-g§) 

+ c(-j(2j - 1 + 2k)q\ + i(2i - 1 + 2k)q 2 2 )) . 

Collecting similar terms in this sum, we find that it van- 
ishes as long as 

aij,k = 5j+ij-i,fe Xj/Xi+i (32) 

for any i, j and k such that i + j + k = n, where we have 
denoted 

Xm = (c(2n + 2 - 3m) + 6(m - l))/m. 

Relation (32) can be regarded as a recurrence for co- 
efficients ciij t k for a fixed k. For k = 0, we start the 
recurrence assuming a„,o,o = 1- This yields 

nn—i y-ri 
m=l Xm H m =i Xm 
u>i.n— 1,0 — TT n 

llm=l Xm 

We obtain now the starting values for the recurrence (32) 
for k > setting, in view of (31), 

O-n-k,0,k = 0-n~k,kfi, 



and find 

~ _ IIro=l Xm Y\jn=\ Xm IIm=l Xm ,_„, 

U l.j,k — yrn ■ \ 00 > 

llm=l Xm 

Evidently, coefficients (33) satisfy the symmetry condi- 
tion (31). Hence, (30) with the coefficients (33) is a sym- 
metric homogeneous polynomial of degree 2n, whose Hes- 
sian commutes with the Hessian of p\ (q). The potential 

$(q,t) = M2(i)^ + E ^ n (t)p^(q) (34) 

n>2 

defines an omni-potcntial flow of a non-Zeldovich type 
in M 3 (provided the coefficients /J,2n{t) tend to zero suf- 
ficiently fast to guarantee convergence of the series (34) 
and to allow its termwise differentiation). Since the poly- 
nomial j4™ (q) is convex for 

< c< 6(n- l)/(n- 2), 

the potential (34) is convex if all /i2n(t) are non- negative 
and < c < 6 (or for < c < 6(N - 1)/(N - 2), if all 
H2n(t) vanish for n > N). 

Although we have constructed our example without 
prescribing the invariants, it might be of interest to cal- 
culate them for the solutions that have been obtained. 
Straightforward calculations yield the values of the in- 
variants for the potential <&(<?, i), for instance, 

( 3,i) = (6 + 3g) g f-(6 + g)^ + 2q 2 (ql-ql) 
721 2Zqiq 2 qi{q% - qfj ' 

which is consistent with a non-trivial dependence of the 
eigendirections on the trajectories (labeled by the La- 
grangian coordinates). We note that, although the solu- 
tion is symmetric in the spatial coordinates, the symme- 
try is lost in the invariant. This stems from the invariant 
under consideration being a nonlinear function of the pro- 
jections of the eigendirections on the plane (^1,92), and 
also from the components of the eigendirections not be- 
ing invariant under all permutations of coordinates (an 
eigenvector is invariant under a permutation of the spa- 
tial variables qi -H> qj provided its ith and jth components 
are swapped). 

We have used two approaches for constructing exam- 
ples of omni-potcntial flow. In three or more dimensions, 
we used the building-block strategy in which the field of 
eigendirections of the commuting Hessians is character- 
ized by prescribing one of the blocks (in our construction 
the polynomial (25)). The problem of commutation of 
Hessians then reduces to three linear equations in the 
unknown block. These equations must be satisfied si- 
multaneously, and we found that in general no solution 
exists for an arbitrary prescribed block. In two dimen- 
sions we have followed another approach, whereby the 
field of eigendirections of the commuting Hessians is char- 
acterized by prescribing the set of invariants, from which 
the field of eigendirections of the commuting Hessians 
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can be uniquely determined. In M 2 just one such invari- 
ant, g(q) (see (12)), should be considered. In R 3 one 
must consider three invariants, for instance, (A7)-(A9) 
of the Appendix, giving rise to three nonlinear equations 
(21)-(23). As in the former approach, these equations 
must be satisfied simultaneously, and hence a solution 
does not exist for an arbitrary set of prescribed invari- 
ants. Thus, whichever approach is used for construction 
of omni-potential flows in K 3 , the prescribed data must 
be tuned for a solution to exist. In the former approach 
the equations are linear and thus simpler, but one is left 
with just one function which can be tuned to achieve 
consistency of the three equations under consideration; 
in the latter approach the equations are nonlinear and 
thus more involved, but one has the freedom of tuning 
three a priori independent scalar functions to gain con- 
sistency of the equations. Of course, in three dimensions, 
the three invariants of the Hessian of the potential can- 
not be prescribed as arbitrary functions. In other words, 
some conditions on the invariants must hold for the three 
equations (21)-(23) to be compatible. Can such condi- 
tions on the invariants be expressed in more explicit form 
remains an open mathematical problem. 



IV. A WKB APPROACH TO 
TWO-DIMENSIONAL OMNI-POTENTIALITY 

So far we have obtained special cases of non-Zcldovich- 
type omni-potential flows. How general are they? Can 
we, for example, in the two-dimensional case prescribe 
an arbitrary smooth initial velocity potential (po(q) or, 
more precisely, the invariant of its Hessian W.(ip (q)): 



linear first-order PDE: 



9f 2 )<A)(g) 



Hn — W.2 



dl2<A)(<?) 



H 



12 



(35) 



that appears in the general equation (12)? In this section, 
where we use only Lagrangian coordinates, d\ and d 2 are 
short for d qi and d qi ; similarly, , d\ 2 and <9f 2 denote 
the second Lagrangian derivatives. 

We shall now show how the construction of non- 
Zcldovich-type omni-potential flow with arbitrary invari- 
ant function g{q) can be done, using an idea of Arnold for 
solving the linear equation which controls the stability of 
solution to the Euler equation [1]. Rather than trying to 
find the most general solution to (12), we construct a spe- 
cial short-wavelength solution through the WKB ansatz 



cD(qr) 



A (q) + -A 1 (q) + ^A 2 (q) + . . . 

K K 



+ CC, 

(36) 

where the wavenumber k is taken very large and where cc 
stands for complex conjugate (needed because we want 
real solutions). In WKB parlance, S(q) is called the 
eikonal function and the functions A (q) 7 Ai(q), ... are 
the amplitudes. 

To the leading order, 0(k 2 ), the WKB ansatz turns 
the linear second-order PDE (12) into the following non- 



(0i S) (diS) 



(37) 



It is easily checked that (37) is equivalent to the state- 
ment that, in the leading order, V q S(q) is an eigenvec- 
tor of the Hessian H(ifo(q)). Actually, this can be seen 
directly, by an argument which also applies in space di- 
mensions d higher than two. Assume that the leading 
WKB term for the potential has a fast spatial depen- 
dence involving the phase factor e 1KS ^ q \ then the Hes- 
sian will involve in the leading order a matrix factor 
—K 2 (diS)(djS). This is a degenerate matrix with one 
eigenvector of non-vanishing eigenvalue in the direction 
of V q S(q); all perpendicular vectors are associated with 
the eigenvalue zero, which has multiplicity d — 1. Omni- 
potentiality requires that this degenerate matrix com- 
mute with the Hessian of the initial potential or, equiva- 
lent^, that VqS(q) be an eigenvector of V.(fo(q))- 

Returning to the two-dimensional case, we now con- 
struct the eikonal function S(q). This construction will 
be done only locally in a neighborhood il, in which the 
Hessian 'H{<fo{q)) is sufficiently smooth and its eigen- 
values are everywhere distinct. (We recall that a double 
eigenvalue is an event of codimension two, which typically 
takes place at isolated locations.) Let n^(q) and n^ 2 \q) 
be two unit eigenvectors of H(ipo(q)), chosen to depend 
smoothly on q in fi. The condition that the gradient of 
the eikonal function be parallel to an eigendirection can 
now be expressed as 



n^>(q)-V q S(q) = or n^\q) ■ V q S{q) = 0. (38) 

In words, these equations state that either n^(q) or 
«.( 2 )(q) is normal to the level lines of the eikonal func- 
tion. Equivalently, the level lines of S are the integral 
curves defined by either n^(q) or n^(q). These form 
a set of orthogonal curves. We thus have two classes of 
solutions. We can prescribe S arbitrarily on one of these 
curves C and extend it locally by demanding that it re- 
mains constant on all the curves orthogonal to C. Note 
that these orthogonal curves play here the role of rays 
in geometrical optics and are thus conveniently called 
"rays" . 

Next wc write the equations for subleading corrections 
obtained by substituting (36) in (12) and identifying the 
coefficients of the various positive and negative powers 
of the large parameter k. We shall only write the equa- 
tions appearing at orders n 1 and n° (the higher-order 
equations have a similar structure). For what follows it 
is convenient to use the compact notation introduced by 
Monge in his theory of surfaces: p, q, f, s and t stand re- 
spectively for d 2 S, d^S, d 2 2 S and d 2 2 S. (We added 
hats to avoid possible confusions.) Furthermore, we write 
g for g(q). The leading-order equation, (37), repeated for 
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convenience, and the two first subleading equations are: 

f-f- gpq = 0, (39) 
(f - gs - i)A + 2(pd 1 A - qd 2 A ) - g(pd 2 A + qd^) 
= 0, (40) 
(f - gs - i)A 1 + 2(pd 1 A 1 - qd^) - 0(^2^1 + qdiAt) 
-i(^i 2 i - dh - gdf 2 )A = 0. (41) 

Using (39) to eliminate the function g, we can rewrite 
(40) and (41) as 

qdl A - P 2 A + ^-'l-Jf-^ A, = 0, (42) 

pz _|_ q 2 

qdxAx - pd 2 A 1 + — 1 — Ai 



p2 _|_ ^2 



i ™ „ 2 (^ii -dl 2 -gd\ 2 



)Ai = 0. 



(43) 



Equation (42) is a first-order linear homogeneous trans- 
port equation for the amplitude y4 along the rays. It 
can be integrated starting from arbitrary non-zero data 
on any curve orthogonal to the rays. Equation (43) for 
the amplitude A\ is of the same sort, except that it has 
an inhomogeneous term involving A . We may thus take 
vanishing data for A\ on an arbitrary curve orthogonal 
to the rays. Higher-order amplitudes satisfy similar in- 
homogeneous transport equations. 

Now we construct locally in space and time an omni- 
potential flow having a given invariant function g(q). We 
take the initial potential (fo{q) arbitrary, but sufficiently 
smooth. Hence, by the WKB method described above 
we can construct a smooth eikonal function S(q) and 
smooth amplitude functions Ao(q), A\(q), .... This, in 
principle, yields a smooth solution, $(q), to (12). (We 
shall not address here the issue of the convergence of the 
WKB series (36).) Because of the imaginary exponen- 
tial dependence on k, the potential Q(q) has very large 
second spatial derivatives 0(k 2 ) and has no reason to be 
convex. However, the following time-dependent potential 
defines an omni-potential Lagrangian map: 

$(<?, t) = ^ + tMQ) + (44) 

Here f(t) is an arbitrary smooth function of time that 
vanishes, together with /, at t = 0. For example, we 
can take f(t) = t 2 . For sufficiently small t and suffi- 
ciently small e, the last two terms in the r.h.s. of (44) 
will not spoil the convexity of the first term. We have 
thus constructed (locally) omni-potential flows for a quite 
arbitrary initial velocity. For large k, the trajectories re- 
sulting from (44) differ only minutely from the straight 
Zeldovich trajectories associated with the two first terms. 
However, these flows are not of Zeldovich type because 
of the third term in the r.h.s. 

If we try to extend the above WKB procedure from two 
to three dimensions, we encounter an obstacle already in 
constructing the eikonal function S. As we have seen, its 



gradient with respect to the Lagrangian position q should 
be everywhere parallel to an eigendirection of the Hessian 
H(ipo)- Denoting by n(q) a unit eigenvector, taken with 
a smooth q-dcpcndcncc, we should then have 



V q S(q)=n(q)n(q), 



(45) 



where p(q) is a scalar function. In other words, the 1- 
form n(q) ■ dq should have an integrating factor (a factor 
which makes it an exact 1-form). This is in general pos- 
sible (locally) in two, but not in three dimensions. 



V. COSMOLOGICAL IMPLICATIONS 

So far our point of view has been kinematical: we 
constructed omni-potential flows without any underly- 
ing dynamical equations. In cosmology the dynamical 
setting is rather well known and discussed for exam- 
ple in Ref. [2]. Let us just recall a few salient points. 
The most widely accepted explanation of the large-scale 
structure seen in galaxy surveys is that it results from 
small primordial fluctuations that grew under gravita- 
tional self-interaction of collisionless cold dark matter 
particles in an expanding universe. The evolution of the 
mostly collisionless matter in the Universe is described by 
the Vlasov-Poisson system in the position- velocity phase 
space. At early times, i.e., close to the epoch of matter- 
radiation decoupling, the expansion of the Universe se- 
lects a single velocity solution at each position rather 
than a distribution of velocities. This feature persists un- 
til particle crossing ("shell-crossing" in the cosmological 
language), where multi-stream solutions are developed. 
As long as multi-streaming is ruled out or is confined 
to scales sufficiently small to be neglected, the Vlasov- 
Poisson system may be replaced by the Euler-Poisson sys- 
tem. Following the notation of Ref. [7] and using Eulerian 
comoving coordinates, x 7 together with a time variable 
r based on the amplitude factor of the growing mode of 
linear theory, we can write the Euler-Poisson system as 



d. 



a 



> T u + (v ■ V x )v 
-p + V x • (pv) 



~27^ + Vx< ^)> 
0, 

P-1 



(46) 
(47) 
(48) 



Here, v is the peculiar velocity, p the density (suit- 
ably normalized) and ip g the gravitational potential. As 
t — ?• 0, to avoid singularities, the density must approach 
unity everywhere; thus the distribution of matter is in the 
leading order uniform as t — > 0. Similarly, v — > — V x <y9 g 
as t — > 0; thus the initial velocity is potential, but other- 
wise arbitrary. It then follows from (46) that the velocity 
stays potential at any later time. 

Reconstruction handles the Euler-Poisson system as 
a two-point boundary-value problem in which the ini- 
tial density is prescribed as uniform and the final (cur- 
rent) density is given by astronomical observations. This 
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is a mass transport problem whose cost function is the 
action associated with the Euler-Poisson equations (see 
Refs. [7] and [17]). Unfortunately, because this action is 
a rather complicated functional, we do not yet possess 
efficient numerical algorithms allowing us to solve this 
mass transport problem. 

The situation simplifies with the Zeldovich approxi- 
mation, which amounts to setting to zero the r.h.s. of 
(46). The remaining equations are then (i) the three- 
dimensional inviscid Burgers equation, which implies 
that particles are moving with constant velocity (in the 
coordinates here chosen) , and (ii) the continuity equation 
expressing mass conservation. The action to be mini- 
mized for reconstruction reduces then to its kinetic en- 
ergy term. As a consequence, the cost function is just 
the mass-weighted integral of the squared displacement 
of fluid particles from their initial (Lagrangian) positions 
q to their current (Eulerian) positions x, as required 
by a theorem of Brcnier [6]. After discretization, this 
mass transport problem becomes an assignment problem, 
which can be solved by efficient algorithms (see Sec. 4 of 
Ref . [7] and Ref . [3] ) . This is the essence of the Mongc 
Ampcrc-Kantorovich (MAK) reconstruction method. 

Mohayaee et al. [18] tested the quality of MAK recon- 
struction by applying it to final states of standard N-body 
simulations, that were performed for various random ini- 
tial conditions of cosmological relevance. The authors of 
[18] noted the unprecedented accuracy of the reconstruc- 
tions down to a few megaparsecs. In particular, they per- 
formed comparisons between three different Lagrangian 
maps: (i) the map based on the N-body integration, (ii) 
the map obtained by applying the MAK procedure to the 
current density field, calculated by the N-body integra- 
tion, (iii) the map obtained by applying the Zeldovich ap- 
proximation, starting from the same initial condition as 
for the N-body simulation. The conclusion of their com- 
parisons is that the N-body map is approximated much 
better by the MAK-generated map than by the Zeldovich 
map. This is particularly striking in their Fig. 7, which 
shows the negative Lagrangian divergence of the displace- 
ment x — q, obtained by the three methods mentioned 
above. 

Can we understand this good performance of MAK re- 
construction on sufficiently large spatial scales? First, 
let us make the rather obvious observation that any 
Lagrangian map that is the (Lagrangian) gradient of 
a convex function (here called the "Brcnier property") 
will be reconstructed exactly (in a discrctized version), 
if we solve the associated quadratic-cost optimal trans- 
port problem, for example, by using the MAK proce- 
dure. Here is a trivial example of this: if we let an 
initially quasi-uniform mass distribution evolve by pure 
Zeldovich dynamics to a final distribution and there is no 
shell crossing, then the MAK reconstruction is exact. If 
the solutions to the Euler-Poisson equations had the Bre- 
nier property, MAK would perform an exact reconstruc- 
tion, but they don't: flows which solve the Euler-Poisson 
equations have no Eulerian vorticity but do generate La- 



grangian vorticity [10]. 

However there exists a refinement of the Zeldovich ap- 
proximation which possesses the Brenier property. This 
is given by the second-order of the Lagrangian perturba- 
tion theory [4, 8, 9, 12, 20, 24]. Here we shall not describe 
the Lagrangian perturbation theory in any technical de- 
tails since the reader can find them in the above pub- 
lications. Nevertheless, in order to discuss some of its 
conceptual problems, we describe briefly a few key steps. 
One rewrites the Euler-Poisson equations in Lagrangian 
coordinates, to obtain a set of nonlinear equations for the 
displacement x — q and its space and time partial deriva- 
tives. Assuming then that in a suitable sense (see below) 
the displacement is small, 0(e), one expands the equa- 
tions in powers of the small parameter e. Here, the only 
perturbed quantities are the deviations of the particle 
trajectory from the homogeneous Hubble flow, i.e., from 
a purely expanding Universe. At the first order, 0(e), 
one has the Zeldovich approximation, which, as discussed 
in Sec. II C, is omni-potential. In particular, the La- 
grangian map is potential (by the Brenier property) and 
the velocity is potential in both Eulerian and Lagrangian 
coordinates. As we have said, with the Lagrangian per- 
turbation theory, one may refine the Zeldovich approxi- 
mation. The second-order Lagrangian perturbation the- 
ory (usually denoted by L2) captures significant gravita- 
tional physics, for example some tidal effects [11], whose 
importance in the large scale structure formation has 
been widely recognized (see, e.g., [23]). L2 has the re- 
markable property that, for standard cosmological initial 
conditions, the Lagrangian map is still potential. As a 
consequence, the velocity is also potential in Lagrangian 
coordinates. With L2, in Eulerian coordinates the veloc- 
ity is potential only up to second order. One would have 
to sum the whole series to arbitrarily high orders of the 
Lagrangian perturbation theory (i.e., to arbitrary orders 
of e) to recover the Eulerian potentiality of the veloc- 
ity (but note that convergence of the asymptotic series 
is not guaranteed [21, 24]). L2 having the Brenier prop- 
erty, the Lagrangian and inverse Lagrangian maps can be 
reconstructed exactly as an optimal transport problem, 
for example, by the MAK technique. This is probably 
the main reason why MAK performs well (at sufficiently 
large scales). Beyond the second order of the Lagrangian 
perturbation theory, scales below the non-linearity scale 
are expected to play a decisive role and Lagrangian vor- 
ticity is unavoidable [9]. Such scales cannot be handled 
accurately by standard MAK reconstruction. 

Finally, let us discuss briefly the thorny issue of the 
validity of the Lagrangian perturbation theory. As men- 
tioned before, the Euler-Poisson equations are a conse- 
quence of the Vlasov-Poisson equations only as long as 
multi-streaming is absent. The problem is that, with any 
cosmologically realistic initial condition at decoupling, 
multi-streaming appears immediately or, anyway, well 
before the present epoch. The situation is somewhat sim- 
ilar to what we would have with a one-dimensional Burg- 
ers flow in which the initial velocity would be spatially 
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non-diffcrcntiable or nearly so, for example, a velocity 
field whose space dependence is the Ornstein-Uhlenbeck 
process: shocks would then appear after an arbitrar- 
ily short time. This, of course, constitutes no reason 
to discard the Euler-Poisson equations for cosmology. 
One way to keep them is to regularize the initial veloc- 
ity field by applying a low-pass filter, i.e., by setting to 
zero all Fourier amplitudes whose wavenumbers exceed 
a cutoff K = l/L, where L is the regularization scale, 
ft may then be shown that no shell-crossing and thus 
no multi-streaming will take place for at least some fi- 
nite time T+(L). Roughly, T±(L) is the inverse of the 
largest shear (velocity gradient) of the regularized initial 
velocity. Up to this time the Euler-Poisson equations 
are valid. Furthermore, it may be shown that, up to 
time T, the small parameter which controls the validity 
of the Lagrangian perturbation theory (including that of 
the Zeldovich approximation) is T/T±(L). Of course, the 
true problem is not regularized, or only barely so. Com- 
parison with N-body simulations suggests nevertheless 
that the Lagrangian perturbation theory remains valid 
for T/T+{L) <C 1 when looking only at scales greater 
than L. Understanding this from a theoretical point of 
view is a challenge, which has been partially addressed 
by Buchert [10]. 



VI. CONCLUDING REMARKS 

The main question that we have addressed in this pa- 
per concerns omni-potentiality, the (convex) potential 
character of the mapping from any time t\ to any time 
t 2 > ti with < ti < t 2 < T. First, we have considered a 
class of flows of "Zeldovich type" , comprised of pure Zel- 
dovich/Burgcrs flows and those obtained from them by 
application of arbitrary nonlinear transformations of the 
time variable and arbitrary time-dependent scale factors. 
Such flows are trivially omni-potential. So are spherically 
symmetric flows. We then have investigated (i) the exis- 
tence of non-trivial omni-potential flows, (ii) their gener- 
icity: can we prescribe the initial velocity potential in an 
arbitrary way? 

The flows have been characterized by their Lagrangian 
maps q M> x{q 1 t) = Vq$(q,t) in terms of the scalar 
potential $. As shown in Sec. II A, omni-potentiality 
implies that along any particle trajectory the Hessians 
H( C E>(<7, t)) commute and thus have common cigcndircc- 
tions. The field of such eigendirections is prescribed as 
a function of the initial (Lagrangian) position q, for ex- 
ample, by the eigendirections of the Hessian of the initial 
velocity potential cf>. The set of eigendirections of a real 
symmetric d x d matrix depends on d(d — l)/2 parame- 
ters and can be characterized, for example, by d(d— l)/2 
of the invariants discussed in the Appendix. As we try 
to determine a single scalar function $, the situation is 
rather different in two and higher dimensions. 

When d — 2, we have a single invariant expressible as 
a ratio of suitable combinations of spatial second deriva- 



tives of $. Thus, prescribing the field of the invariant 
values, g(q), we obtain a linear second-order PDE, (12), 
for $. For a suitable family of fields g(q), we have found 
in Sec. Ill A non-trivial omni-potential flows that are lin- 
ear combinations of homogeneous polynomials, thus en- 
suring the existence of such non-trivial flows. Using a 
WKB method, in Sec. IV we have then been able to con- 
struct omni-potential flows, at least locally in space-time, 
for arbitrary smooth g(q). These flows are actually close 
to Zeldovich flows with straight trajectories (but are not 
of "Zeldovich type"). Extending this construction glob- 
ally in space and avoiding the rapid spatial oscillations 
inherent to a WKB method constitute interesting open 
problems. 

When d = 3, omni-potentiality can be expressed in ei- 
ther of two equivalent ways. One is to demand the com- 
mutation of the Hessians of $ with those of a prescribed 
4>; this gives d(d— l)/2 linear homogeneous second-order 
PDEs. The other involves working with the invariants, 
introduced in the Appendix, which are rational functions 
of the entries of the Hessians of $; this gives d(d — l)/2 
nonlinear second-order PDEs. In both approaches we 
have one unknown scalar function, $, which has to sat- 
isfy more than one equation. Hence, there is an issue of 
compatibility of these equations. However, by restrict- 
ing the potential $ to possess a suitable finite symmetry 
group, we have obtained a fairly large class of non-trivial 
solutions that are even-degree homogeneous polynomials. 
Whether non-Zeldovich-type omni-potential flows exist 
for an arbitrary smooth <j> remains an open problem. In 
dimensions d > 3 the situation is basically the same. 

We have shown in Sec. II B that omni-potentiality of a 
flow is equivalent to having at each time a velocity field 
that is potential in both Eulerian and Lagrangian co- 
ordinates. Such double potentiality was frequently con- 
sidered in cosmology. It is of particular relevance when 
performing reconstruction by convex optimization tech- 
niques such as the Monge-Ampere-Kantorovich (MAK) 
procedure. Note that the Euler-Poisson flow is poten- 
tial in Eulerian coordinates, but not in general in La- 
grangian coordinates. As discussed in Sec. V, the ap- 
proximate Euler-Poisson flow obtained by the second- 
order Lagrangian perturbation theory (L2) is exactly po- 
tential in Lagrangian coordinates — and thus its inverse 
Lagrangian map can be obtained exactly by the Monge- 
Ampere-Kantorovich (MAK) procedure — however, it 
is only approximately potential in Eulerian coordinates 
and thus must be qualified as an approximately omni- 
potential flow. In particular, it does not represent an 
example of an exactly omni-potential three-dimensional 
flow for an arbitrary smooth initial velocity. 

We finally wish to mention a concrete open problem of 
cosmological interest: as mentioned, MAK gives the ex- 
act inverse Lagrangian map for L2 (and this contributes 
to explaining why MAK works so well when tested with 
N-body simulations). However, in L2, between the Eule- 
rian position x and its Lagrangian antecedent q, the tra- 
jectory is not given exactly by the Zeldovich approxima- 
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tion that has a constant velocity (x — q)/r (in the coor- 
dinates we used in Sec. V). The actual L2 trajectory is 
slightly curved and its current (peculiar) velocity differs 
slightly from (x—q)/r. It would be of interest to find how 
to perturbatively handle such discrepancies. Given that 
the full Euler-Poisson reconstruction problem and the 
Zeldovich approximation to it both have convex optimi- 
zation formulations, the question arises, whether L2 and 
higher-order approximations possess such formulations. 
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Appendix: Invariants under variation of eigenvalues 
of symmetric matrices 

Here we show how the set of commuting symmetric 
real d-dimensional matrices, having prescribed eigendi- 
rcctions, can be characterized by a certain number of 
invariants, which are rational functions of the matrix en- 
tries. The findings here may be of interest beyond the 
study of omni-potcntial flow. To the best of our knowl- 
edge, these results are not available in the published lit- 
erature. If the reader is aware of any relevant reference, 
kindly inform the authors. 

By the theorem on codiagonalizability of commuting 
symmetric matrices, in an omni-potential flow, Hessians 
of the potential of the Lagrangian map must have the 
same set of eigendirections at all times along any trajec- 
tory. We therefore ask ourselves a more general question: 
Suppose, the eigenvalues of a symmetric matrix are be- 
ing varied while the eigendirections remain fixed. Which 
quantities constructed from the entries of the matrix are 
unaltered under such variations? We call such quanti- 
ties invariants. For a two-dimensional matrix Hij , as we 
have seen in Sec. II, the eigendirections depend on the 
single invariant g — (Hn — Hyi) I When d < 5, our 
question can, in principle, be answered by first solving 
the characteristic equation whose roots are the eigenval- 
ues and then determining the eigendirections. However, 
the expressions for the eigendirections will then involve 
radicals, so that obtaining rational expressions for the 
invariants is not easy. Anyway, when d > 5, the non- 
solvability of the characteristic equation by radicals ren- 
ders this strategy useless. We need therefore a more prac- 



tical general algebraic approach to this problem. 

1. Construction of invariants in dimension d 

Let us assume that all eigenvalues Aj of a symmetric 
d x d matrix H are distinct, and thus all eigendirections 
are uniquely defined. Description of an arbitrary set of 
d orthogonal eigendirections in R d requires d(d — l)/2 
parameters: d arbitrary directions require d(d — 1) pa- 
rameters, from which one must subtract the number of 
orthogonality conditions, d(d — l)/2. We expect there- 
fore that a set of d(d — l)/2 suitably chosen invariants 
uniquely defines the eigendirections. 

The problem of finding such invariants is an instance 
of a much more general problem of characterizing lin- 
ear subspaces of a vector space; here the vector space 
is that of all real symmetric matrices and the subspace 
that of matrices having prescribed eigendirections, which 
is spanned by the set of all the powers, from zero to d— 1, 
of this matrix. The general problem can be, in principle, 
handled using Pliicker coordinates [14]. For our problem, 
a more direct approach is available, as now explained. 
For d > 3, our characterization involves fewer invariants 
than the corresponding number of Pliicker coordinates. 

We denote by h(Aj) an eigenvector associated with the 
eigenvalue Xi , and assume without any loss of generality 
that no component of any eigenvector vanishes; one can 
always achieve this by suitably rotating the orthonormal 
basis in R d , in which the eigenvectors are decomposed. 

We construct the invariants as follows: We set, for 
some l<m^n<d and k < d, 

7 { ± k) =P id ' k) Wmn,l,...,Pmn,d), (Al) 

where 

flmn,i = h m (Xi) I h n (Xi) 

and p( d,fc ) denote symmetric homogeneous polynomials 
of degree k < d, 

p(«W( y ) = Yl Vh-Vii-Vi* 

i<h<-<n<-<ik<d 

for y e R d . By construction, the quantities 7™' de- 
pend only on the eigendirections (through the ratios of 
components) and are invariant under permutations of 
the eigendirections; thus they depend only on the set of 
eigendirections. Then one substitutes into (Al) the re- 
spective components of the eigenvectors h(Aj), expressed 
in terms of the associated eigenvalues A, and of the entries 
of the matrix H. It is easily seen that this will produce ra- 
tional functions of the matrix entries and of the eigenval- 
ues. Furthermore, it may be shown that the eigenvalues 
enter only through symmetric polynomial combinations, 
which — by Viete's theorem applied to the characteristic 
polynomial — have a polynomial dependence on the ma- 
trix entries. The actual derivation of the invariants can 
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be partially simplified by making use of the identity 



H(\ k + c) = det\\H + ci\\. 



(A2) 



fe=i 



2. Relations between invariants 

We obtain thus d 2 (d— 1) invariants ~fmn^ in the form of 
rational functions of the entries of the symmetric matrix 
H. Evidently, these invariants are too numerous to be 
independent. For instance, for any l<m^n^l<d 
and < k < d they clearly satisfy the relations (no 
summation on repeated indices!) 



(d,d) (d,d) 

Imn inm 



1. 



and 



(d,d) (d,d) 

'ml 1 In Imn 



Ad,k) = ld,d) ld,d-k) 
Imn Imn inm 



(A3) 



(A4) 



(A5) 



Identities (A3) and (A5) link invariants "fmn for, say, 
m < n with those for m > n; there are d 2 (d — l)/2 
of such independent relations between invariants. Equa- 
tions (A4) imply 



m—l 

(d,d) _ IT Jd,d) 



11 7»+l,i 



for any m > n + 1; conversely, any relation (A4) follows 
from these relations together with (A3). Thus, the iden- 
tities (A4) contribute further (d— l)(d— 2)/2 independent 
relations. 

For any n such that 1 < n < d, the relation 



imply, for each p > 0, the relation 



d d 



/ j / j mn Imn 



n—1 m—l 

m^n 

d d d 

V V V hW— 

2^ 2-* ,2-, mn h (A 

t—l n—1 m—l v 

m^n 

EE(Af-^) = o. 

i—l n—1 



(A6) 



Since, by the Cayley-Hamilton theorem, any matrix 
is a root of its characteristic polynomial, the entries 

H&1 for p > d are linear combinations of Hmn for 
p — d < p' < p — 1, the coefficients in these linear combi- 
nations being independent of indices p, m and n. There- 
fore, relation (A6) for p > d is a consequence of d such 
relations for p = p' such that p — d < p' < p — 1. The 
relation (A6) for p — is trivial, and hence there are d— 1 
independent relations ( A6) for 1 < p < d — 1 . 

In principle, the total number d 2 {d — 1) of invariants 

Jmn should exceed the number d{d — l)/2 of indepen- 
dent invariants by the number of relations between the 
invariants. For arbitrary d, we have obtained above 



d 2 (d 



1) , (<* -!)(<*■ 



2) 



+ d+(d-l) 



rf(d 2 + 1) 



independent relations. For d = 3, they constitute 15 re- 
lations constraining the 18 invariants that we have intro- 
duced; this fits our expectations that at most 3 invari- 
ants are independent, because the orthogonal frame of 
eigendirections of a 3 x 3 symmetric matrix is described 
by 3 Euler angles. For d > 3 the number of the obtained 
independent relations is still insufficient to fill the gap, 
and more relations are to be identified. 



E 

m—l 

m^n 



imn 



d(d- 1) 



stems from orthogonality of the eigendirections (there are 
d such relations). 

Another family of relations involves the invariants 
7™' '. For p > 0, let denote the pth power of the 
matrix H and let Hmn denote its entries, which are of 
course readily expressed in terms of the entries of the 
matrix H; we also set = I. The relations 



m—l 



that hold true, for any i and n, by definition of eigenvec- 
tors of H, h(Xi), and the identity 



EA? =trHM 



Invariants for d = 2 



As an example, we present a detailed derivation of the 
(2 1) 

invariant 7^' ' for d = 2. The ith eigenvector of a 2 x 2 
symmetric matrix H is (i?i2, \ — Hn), and hence 



(2,1) 
Tl2 



H 



12 



+ 



H 



12 



H12 (Ai + A 2 — 2i?n) 
Ai — iiii A2 — -ffll (Ai — iin)(A2 — -ffn) 

The characteristic equation for the eigenvalues is 

A 2 — (Hn + i?22)A + H11H22 — H 2 2 = 0, 

and hence Ai + A 2 = H\\ + H 2 2- By virtue of (A2), 



(Ai- J ffn)(A 2 -ff 11 ) = det 
Consequently, 



H 12 
H12 H22 — Hn 



- ~Hi2- 



(2,1) Hu — H22 
I12 = 77 • 

-"12 
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The invariant 7^' ^ can be found by a similar calcula- 
tion, or just by swapping subscripts in the expression 
for 7 1 2' 1 ' 1 ; clearly, the two invariants are interrelated: 

721'^ = — 7i2' 1 ' ) - The invariant 721'^ turns out to be 
degenerate: 



(2,2) 
721 



H 12 



(Ai - JTn)(A 2 - Hu) 



= -1. 



Thus, we have obtained the same invariant as that found 
in Sec. II D. 



4. Invariants for d — 3 

We consider now invariants in the three-dimensional 
space. The zth eigenvector of a 3 x 3 symmetric matrix 
H has components 

(Hi 2 H 2 3 + Hi 3 (\ - H22), Hi 2 Hi 3 + H 23 {\ — Hu), 

(Ai — H 1 i)(X i — H22) — H 12 ). 
The procedure outlined above yields 
(3,1) H22 — Hu 

721 = 77 

-"12 



+ 



H13 (Hu — H 22 )H\ 3 H 23 + {H 23 — H 1Z )H\2 
H\2 (H22 — H 33 )H\ 2 H\ 3 + (Hf 3 — Hi 2 )H23 



[Hu — H 33 )Hi 2 H 23 + {H 23 — H 12 )Hi3 

(H22 — H 33 )Hi 2 H\ 3 + (H 13 — H 12 )H23 



(A7) 



(3 3) 

7 21 ' is the ratio of two polynomials, which we calculate 
using (A2): 



(3,3) _ (#11 — H 33 )Hi 2 H 23 + (H 23 - H^ 2 )H 13 
(H 22 — H 33 )Hi 2 Hi 3 + (H 13 - H 12 )H 23 



721 = 



72i' 2 ^ can be found by applying identity (A5): 



(3,2) (3,3) (3,1) 
721 = 721 7l2 



(A8) 



(A9) 



(3 1) 

(here 7 12 ' can be obtained by permuting the subscripts 
1 and 2 in (A7)). 

The three invariants 7 2 i' fe ^ for 1 < k < 3 uniquely 
define the three ratios 02i,i'- by Viete's theorem, they 
are roots of the cubic equation 

/3 3 - 7 r ) /3 2 +7^V-7r ) =0. 

The eigendirections can be recovered in the form of three 
eigenvectors (1, /5 2 i,i, q). One can try to obtain the third 
components q (i = 1, 2, 3) from the relations express- 
ing orthogonality of the eigendirections. However, this 
produces two solutions: if {q} is obtained in this way, 
then {— Ci} is also a solution. Hence, the invariants 721'^' 
1 < k < 3, define two distinct sets of eigendirections. The 



non-uniqueness is eliminated, if in addition we know any 
of the invariants 7^3'^ or 7jy'^ for i = 1,3 and j = 1,2. 
It is unclear, whether one can choose a set of three in- 
variants uniquely defining three eigendirections. 

When all invariants "fmii ^ are known for k = 1 and 3, 
the equations for the entries of symmetric matrix H can 
be considerably simplified. In view of (A8) and the same 
equation with permuted subscripts 2 and 3, relation (A7) 
can be expressed as 



H 



22 



-"12 -"12 



Adding (A10) to its analogue, where subscripts 1 and 2 
are permuted, we obtain 

-H a (7^ 1) +7r ) +7& 1) +7f 2 3) ) 

' "l3-:;f • "2::-;;f ». (All) 

Permuting subscripts in this equation, we obtain a lin- 
ear system of equations for the non-diagonal entries of 
H. (Here, we note that the sum of (All) and its coun- 
terparts with permuted subscripts reduces to (A6) for 
p = 1.) Upon solving this linear system, we find the 
differences H mm — H nn from (A10) and the analogues 
of this equation with permuted subscripts, i.e., we de- 
termine all diagonal entries up to an additive constant. 
Further determination of the matrix H would require, of 
course, the knowledge of its three eigenvalues. This im- 
plies that the non-diagonal entries, as determined from 
the above-mentioned linear system, involve two free pa- 
rameters. This, in turn, requires that (All) be equiva- 
lent to any equation, obtained from it by permutation of 
subscripts. 

Consequently, we obtain further relations between the 
invariants: Permuting subscripts, say, 1 and 3 in (All), 
we get 



-H 



23 



(3,3) . (3,1) 
' 723 + 732 



^13713 



H\2l\2 



0. 



This equation is equivalent to (All) if and only if 

/ (3,1) . (3,3) . (3.1) . (3.3)\ (3.3) . (3.3) (3.3) n 
(721 + 721 + 7l2 + 7l2 ') 7l3 + 731 7l2 = 

and 



/ (3.1) . (3,3) . (3,1) . (3,3)\ 
(721 +721 +7l2 +7i2 ) 

„ ( (3.1) . (3,3) . (3.1) . (3,3) 
X (723 + 723 + 732 + 732 



\ (3,3) (3,3) 

) = 732 7l2 • 



In view of (A4), the first of these relations is equivalent 
to 

(3.1) , (3.3) . (3,1) . (3,3) . (3.3) (3,3) n / » 1 x 

7 2 i' + 721 + 7i 2 + 7i 2 + 7 3 i 7 32 = 0, (A12) 

and the second one follows from (A12) and the relation 
obtained from (A12) by permuting subscripts 1 and 3. 
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Relation (A12) and its analogues with permuted sub- consistency of the invariants, instead of using (A6), one 
scripts can, of course, be established directly. Such re- of our basic relations between invariants, 
lations can be used in three dimensions for verifying the 
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